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This paper reports on a new algorithm to compute the asymptotic solutions of a linear differential system. 
A feature of the algorithm is the ability to accommodate periodic coefficients. 
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1. Introduction 

In this paper we discuss a new algorithm for estimating and improving error terms in 
the asymptotic solution of linear differential systems. We consider systems of the form 

where Z is an n— component vector, p is a real or complex scalar factor, D is a constant 
q ■ n x n diagonal matrix, 

D = dg(d 1 ,...,d n ) (1.2) 

with distinct d/., and i? is also annxn matrix whose entries tend to zero as x — > oo, that 
is, R(x) — ► as a; — > oo. 
If it is the case that p(x)R(x) is L(a, oo), the Levinson asymptotic theorem flEastham 



X 



1989| , section 1.3);(|L evinsonj |1948|) states that there are solutions Z^ (1 < k < n) of ( |1 . 1|) 



such that 

= {e fc + %(a;)}exp(4 / p(t)dt) : (1.3) 

J a 

where is the unit coordinate vector in the k— direction and %(x) — > as x — > oo. 
The size of the error term r]k is related to the size of R{x) as x — > oo, and therefore the 
accuracy of ( |1.3|) can be improved if the perturbation R(x) can be improved - that is, 
made smaller in magnitude — as x — > oo. Under suitable conditions on p and R, this 
improvement can be effected by applying a sequence of transformations to the solution 
vector Z in ( |1 . 1|) . Our algorithm is concerned with the implementation of this sequence 
of transformations. 
In order to introduce the ideas involved, we consider the transformation 

Z = {I + P)W, (1.4) 
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where / is the n x n identity matrix, dgP = 0, and the off-diagonal entries of P are defined 
by 

PD - DP = R- dgR. (1.5) 
Thus, in terms of the (i, j) entries in the matrices, 

Pij = r ij/ ( d j ~ d i) (« ^ j)- (!-6) 
On substituting ( Ofi into ( |1.1[ ) and using ( |1.5| ), we have 

W = p{Z) + (I + P)-\RP - PdgR - p^P'^W, (1.7) 

where 

D = D + dgR. (1.8) 

By ( |1.6| ), P(x) — > as x — > oo and therefore it is clear that there are circumstances, to 
be detailed later, in which the perturbation term in the I^-system has a smaller order 
of magnitude for large x than the original perturbation R. Repetition of the process 
successively improves the perturbation term. It is to the final system in the process that 
the Levinson theorem ( |1.3|) is applied, when a prescribed accuracy in the error term has 
been achieved. 

The transformation back from the final system to the original system (|1.1|) yields an 
improvement of (|1.3|) in which the factor + %(#) is replaced by 

{I + P (x)}{e k + r] k (x)} (1.9) 

with a new r\ k which has the prescribed accuracy, and the matrix P Q is generated explicitly 
by our algorithm. The terms in P tend to zero as x — > oo and r/k = o(Po). Thus ( |1.9| ) 
provides explicit sub-dominant terms for the asymptotic solution of ( |1.1| ). In sections 2 
and 3, we discuss the sequence of transformations and, in sections 4 and 5, we discuss the 
algorithm for the generation of the terms in Pq. 

In a recent paper ( [Brown, Eastham, Evans fc Mc(Jormack||1996| ), we consider a particular 
system of the form ( |1 . 1| ) which arises from the n— th order differential equation 



y 



^(x)-Q(x)y(x)=0, 



;i.io) 

and we formulated an algorithm which implements a sequence of transformations of the 
The main emphasis in ( |Brown et al.| |1996|) , however, is on the analytic 

and for 



type (OD-dH 



1.10) 



and asymptotic implications of the transformations for the solution of 
applications to spectral theory. Here, on the other hand, we wish to develop our algorithm 
from the point of view of symbolic algebra in the context of the general system ( |1 . 1|) . We 
also demonstrate the versatility of our procedure by applying it to other situations than 
the one covered in ( [Brown et al.| |1996| ) . 

Finally in this introduction, we note that the origins of the transformation ( |1.5| )-( [1~6| ) 
lie in the work of Harris and Lutz ( [Harris fc Lutz 1974 ) with more recent developments 
of these ideas by Eastham ( [Eastham] 1986) ( [Eastham 1989j , section 1.7). The nature of 



the matrix / + P is that it is an explicit approximation to the matrix whose columns are 
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eigenvectors of D + R, these eigenvectors in general only being explicit in terms of D and 
R when n — 2. 



2. The sequence of transformations 

We now define a sequence of transformations 

Z m = (I + P m )Z m+1 (m = l,2,...) (2.1) 
of the type introduced in (|1.4|) - ( |1.6|) , the purpose of which is to produce differential 



systems for the Z m , similar to ( 171]) , but with the perturbation term successively improved. 
The definition is almost the same as that given in (|Brown et al.| |1996| , section 4) for the 
particular system (|1 . 1| ) which arises from ( |1.10|) , and so we shall be brief in this part of 
the paper. A typical system in the process is 

Z' m = p(D m + R m )Z m , (2.2) 

where D m is diagonal, with (|1 . 1| ) being the case m = 1. The process ends at m = M when 
the perturbation R M has a pre-assigned accuracy in terms of its order of magnitude as 
X — ► oo. 

As already indicated by (|1.7|) , R m will contain terms of different orders of magnitude as 
x — > oo, and it is the essence of our algorithm to identify and collate these terms according 
to their size. Hence we write 

Rm = V m + E m = Vim + V-Zm + ••• + V^m + E m , (2-3) 

where 

V km = o(V jm ) (x -> oo, k > j) (2.4) 

and 

E m = o(V^ m ) (x oo). (2.5) 

Here E m represents terms which are already of the pre-assigned accuracy, and they take 
little part in the transformation process fl2.1|) . The Vj m represents terms which are not of 
that accuracy, and they are successively replaced by smaller-order terms as we go through 
the process. Also as indicated by (|1.8|), we take any diagonal terms in Vi m over to D m in 



( |2.2j ). Thus we arrange that 

dgVi m = (2.6) 

and we write 

D m = D + A m . (2.7) 

We substitute ( |2.1| ) into ( ^.2[ ) to eliminate the dominant term V\ m in ( ^.3[ ) and to define 
the resulting terms Vj >m+ i constructively in terms of the Vj m - Corresponding to (|1.5|), we 
define P m by 

P m D - DP m = Vim (2.8) 
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with dgP m = 0. Then it is easily checked that ( |2.1| ) and ( |2.2| ) give 
7* 

J m+l 



Z' m+1 = p{D m +(I + P m )-\ - p- x P' m + T m + V lm P n 



+ (R m -V lm )(I + P m ))}Z m+1 , (2.9) 

where 

T m A m P m P m /\ m . (2.10) 
As in (|Brown et al.| |1996|, section 4) , we show that ( !T9 ) can be expressed as 



Z' m+ i — p(D rn+ i + R m+ i)Z rn+ i (2-11) 

where R m +i has the form 

,m+l 

+ E m+ i (2.12) 

as in (|2.3| )- (|2.5|) , but with a different \i. To do this, we let U denote any of the terms on 
which (J + P m ) _1 acts in (|2.9|) . Then we write 

(/ + Pm)- 1 = I ~ P m + Pi ~ - + (-1)^ + (-l)^^/ + P m )- l P» m +1 (2.13) 

where, for each U, v is chosen so that the product 

{I + P m )- l P v m +l U (2.14) 

has a sufficiently small order of magnitude to be included with E m and form part of E m+ i. 
Now we group together terms of the same order of magnitude and denote the dominant 
term by S m+ i. We then obtain (|2.12|) (with S m+ i in place of Vi )Tri+ i), where E m+1 has the 
pre-assigned accuracy and, by (|2.8|) , S m+ i and the V^ m+ \ are known explicitly in terms 
of the Vj m . Then, finally, we obtain (|2.11|) and . 1 2|) by defining 



D-m+i — 

and 

Vx t m+i — S m+ i — dgS'm+i. 
3. Orders of magnitude 



The transformation process (2.1) is carried out for m = 1,2, ...,M — 1 and, in order 
to express the process in terms of an algorithm which can be implemented in the sym- 
bolic algebra system Mathematica, it is necessary to specify more precisely the orders of 
magnitude involved. The starting point is m — 1 and, in (|2.3|) , we suppose that 

V l = Vn + V 21 + ... + V N1 , (3.1) 

where 

Vjxix) =O(x- i) (l<j<N) (3.2) 
as x —>■ 00 and, corresponding to ( |2.4| ), 

< 61 < 6 2 < ... < 9 N . 
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We assume that the 9j in ( |3.2| ) are chosen to have their minimum possible values and, in 
practice, fl3.2| ) represents the exact order of magnitude of Vj\. We denote by a the set of 
positive numbers 

a = {n^i + n 2 9 2 + ... + n N 8 N ; n x > 1, n 2 > 0, n N > 0} (3.3) 

the rij being integers. It is possible that different values of the rij give the same number in a 
and, allowing for this, we denote the distinct numbers in a by eri, a 2 , ... in increasing order. 
Let us suppose that the pre-assigned accuracy represented by E m in ( |2.3| ) is expressed as 

E m (x) = 0(x~ K ) (3.4) 

for a given K > 0. Then we choose the integer L so that 

ol < K < a L+1 . (3.5) 

The definition of a in ( |3.3j ) allows us to postulate orders of magnitude 

V jm (x) =0(x-^->) (3.6) 

where we can allow the possibility that some of the Vj m ( even Vi m ) may be identically 
zero. To justify (|3.6j), we note first that P m = 0(x~ am ) by ([2.8|). Then, recalling the use 
of ( |2 . 1 3| ) in ( |2^9| ), we also have 

PZ l V jm = 0(x-™™-^-i), 



and again ro m + u m+ j_i G a by ( |3.3|) . Further, since the combination r = and j = 1 
does not occur together here, we have 

The term T m in ( p.9|) is treated similarly. A simple induction argument on m now estab- 
lishes ( |3.6|) for all j and m, provided only that we add a suitable hypothesis on the term 
p~ x P' m which appears in (|2.9| ) but is not so far included in the argument. We therefore 
add the hypothesis that 

p- 1 P' m = W lm + ... + W lm (3.7) 

where, similarly to (|3.6|), 

W jm {x)=0{x-^) (3.8) 

and again we allow the possibility that some Wj m may be zero. Since P' m depends on 
V[ m ( see (|2.8|) ), which in turn depends on the previous matrices in the process ( |2.i|) , the 
nature of (|3.7|) and ( |3.8| ) is that they are conditions on the successive derivatives of the 
original Vj± which occur in Ri in ( |2.3| ) and ( p.l| ). The exact form of these conditions on Vji 
determines classes of matrices R\ to which this theory and the consequent algorithms are 
applicable. Examples of such classes will be given in section 5. Thus ( |3.7| ) and ( |3~B| ) are 
consequences of these conditions on Ri which must be established ( usually by induction) 
in each application of the theory. It is these Wj m which will appear in our algorithms. 
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We can now summarise this section by saying that, subject to fl3.2|) , ( |3.7| ) and 
we have established that 

V jm (x) = Oix-^- 1 ) 

in ( p.3|) -( p3|) . Also, allowing for the fact that some Vj m in (2.5) may be zero, we can write 



fi = L — m + 1 by ( p.4| ) and (|3.5|) . The transformation process (|2.1|) ends when (|2.3|) 
reduces to 

i? M = E M = 0(x~ K ), (3.9) 



the pre-assigned accuracy, and it follows from ( |3.5| ) that 

M = L+l,/j, = M-m. 

4. The algorithm 

In this section of the paper we show how the theory that has been developed in sections 
1 through 3 may be used to obtain a computer code to calculate the asymptotic expansion 



of the solutions of (LI) together with an explicit error bound at some point x > X > 0. A 
consequence of the theory that we have exhibited is that, given sufficient computational 
power, the quality of the asymptotics that we obtain allows us to take X to be quite small 
and still maintain high accuracy in the solutions. 
As in the discussion in ( [Brown et al.||1996|) the algorithm is formulated and implemented 



in three stages. All the symbolic algorithms that we shall discuss are implemented in the 
symbolic algebra system Mathematica. The first algorithm, which is concerned with the 
generation of a set of recurrence relations to compute the matrix quantities Sj, assumes 
only that the quantities involved satisfy non- commutative multiplication. We recall the 
comments made after (|3.8|) that general classes of matricies Ri to which the algorithm is 
applicable will be given in section 5. In the following, we write A m — ( I + P m ) _1 and we 
note that expressions such as P m , T m and Wj m appear in the algorithm by virtue of their 
orders of magnitude as indicated in section 3. 



Algorithm 4.1. ( a ) Define K to specify the accuracy p^j ). 

( b ) Define N and 8i,...,8n in (\3-3j) and arrange the distinct numbers in the set o in 
increasing order. This defines o~i,o~2, ... and determines L in ( \3.d[) . For a given K, nj in 
( P^ ) satisfies 

(I)0< nj <[Z=*-] (j>2) 

(n)i< ni <\^}. 

( c ) Start with D x and V jX (1 < j < N) as in (g^-flE^ and put E 1 = 0. 



( d ) For m — 1 to M — 1, 
( I ) E m+ i = A m E m (I + P m ). 

( II ) For each U G {W jm (1 < j < l),T m , V lm P m , V jm , V jm P m (2 < j < M - m)}, 
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( i ) In l\2. T\) determine v. 



( ii ) For r = to is, 

determine the order a m+ k = ra m + (order of U) of P^JJ ; 
Update V Km+1 = V k>m+1 + {-l) r P m U. 

( in ) Update E m+1 = E m+1 + {-If A m P m +1 U . 

( III ) Output S m+ i = V 1>m+1 . 

At each stage of the algorithm, S m+ i depends on the terms Dk,Pk, Wjk and V k (1 < k < 
m). However, because of part (6), the algorithm requires more precise information than 
its counterpart, Algorithm 6.1, in ( [Brown et al.| |H?9~6"D . The set a in ( jBrown et al.| |T9"9~6] ) 



has a very simple form, consisting only of numbers na, where n is a positive integer and 
a (> 0) is a parameter. Thus a n = na in (|3.3|), and Algorithm 6.1 in ([Brown et al.| [xT^ED 



can be executed without specifying the value of a. We give a more general example of the 
same situation in Example 5.1 below. However, in the wider context of ( |3.3|) , sufficient 
information about the parameters 8±, On must be provided to Algorithm 4.1 in order to 
generate all the necessary values of a n . We therefore defer further discussion of the output 
of Algorithm 4.1 to Example 5.2 in the next section, where values of the parameters are 
specified. 

Algorithm 4.2. Starting with the precise form of the matrices D 1 and V\ and with 
Ei = 0, the expressions S 2 , Sm-i generated by Algorithm are evaluated in order. 
These are then used to evaluate the matrices D m+ i and V\ )m+ i. 



The structure of the algorithm is similar to that of Algorithm 6.2 of ( Brown et al. 1996 ). As 



noted in that paper, in order to reduce the computation time a detailed assessment of the 
mathematical issues involved at each simplification of an expression must first be made. 
Thus the judicious use of the Together, Apart commands instead of the Simplify command 
can result in a dramatic decrease in the time needed to perform the computation. At this 
stage it is necessary to keep the expressions in symbolic form since the Wj m must be 
obtained explicitly. These are computed in terms of P m , which in turn is obtained from 
Sj by differentiation. 

The final algorithm that is needed in the symbolic part of the computation is concerned 
with obtaining an upper bound for the norm || E m || of the error term E m . The norm is 
computed using the sup. norm by applying the triangle and Cauchy inequalities 

|| AB ||< n || A (HI £ ||, || A + B ||<|| A II + II B || 



for n x n matrices. As in ( [Brown et al.| |1996| ) a bound for the norm of the inverse matrix 



A m = (I + P rn ) 1 is given by 

II A m \\< 1+ || P m || /(l-n|| P m ||) (4.1) 

provided || P m ||< 1/n. 
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Algorithm 4.3. Compute the sup. norm of each matrix in E m , using for the in- 
verse matrices. Next apply the triangle and Cauchy inequalities to obtain an upper bound 
for the sup. norm of E m itself. 

We note that, since Algorithm 4.1 expresses E m in terms of matrices arising at earlier 
stages of ( |2.2| )- and therefore ultimately in terms of the Vji in( |l.l| )- so also Algorithm 4.3 



ultimately expresses the norm of E m in terms of norms derived from the original system 

Q)- 

Before moving on to examples of the implementation of the algorithms, we add some 
detail to ( |1.9| ) concerning the generation of the sub-dominant terms in the asymptotic 
solution of (|1.1|). By ( ^.9| ) and (|3.5|), the final system in the sequence ( |2.2| ) is 

Z' M = p{D M + E M )Z M , (4.2) 

where 

II E M \\<c M x- aM (4.3) 

and Cm is a constant. We recall that the numbers a m cover all orders of magnitude which 
occur. Algorithm 4.3 provides a definite value for cm in any particular example. As in ( 
OP , the asymptotic solution of fl4.2|) has the form 



{e k + ri k (x)} exp( d kM (t)p(t)dt), (4.4) 



where the d k M are the diagonal entries in Dm and the size of r\ k can be expressed in terms 
of cm and &m as in (|Brown et al.| |1996| , (3.15)). What we wish to emphasise here is the 



transformation back from ( f4.2|) to the original system (|1 . 1| ) . As indicated in ( |1.9|) , this 
adds to flO] ) the extra factor 

/ + P (x) = U^{I + P m (x)}. (4.5) 



Now the definition of P m in ( |2.8|) is in terms of Vi m , which is provided by Algorithms 4.1 
and 4.2. Further, by and (U), P m (x) = 0{x~ am ) (1 < m < M - 1). Thus, in terms 
of (^4.5|) , our algorithms provide sub-dominant terms up to 0(x~ GM ~ 1 ) in the asymptotic 
solution of (|I~T|) . 

We mention one further point concerning the transformation process which leads from 



(|2.2j ) to (|2.9|) . Since the derivative P' m appears in ( p.9|) and since P m ultimately depends on 
V\ and p, each step in the process requires the existence of a further derivative of V\ and 
p. Thus the sub-dominant terms in ( [4.5|) require the existence of M — 1 derivatives of V\ 
and p. If V\ and p are infinitely differentiable then, subject to convergence considerations, 
I) would yield a full asymptotic expansion. It is hoped to deal with this matter in a 



future paper. 



5. Examples 
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5.1. Example 1 

Let p{x) = x 1 and R(x) = x~( 1+ ^C, where 7 > and C is a constant matrix. Here we 
have just N = 1 in ( |3.1| ) and 

a m = 771(1 + 7) (m = 1,2,...). 

This example is basically the case considered in ( [Brown et al.| |1996| , section 5) with a 
special choice of C and, as in (|Brown et al.| |IT)96l ), the condition ( |3.7|) is easily verified by 
induction on m. The present code has been tested on this example and the results from 



Algorithm 4.1 are, up to notational differences, identical with those reported on in ([Brown 



et al.||I~99~6|) . Further, Algorithms 4.2 and 4.3 return values of the solutions computed with 
4 iterations that, at X = 40, are within 10 -11 of those reported on in ( Brown et al. 1996 ). 

5.2. Example 2 

A significantly different example is obtained when p(x) and R(x) in (|1.1|) contain peri- 
odic factors. Let p(x) = x 7 0(x /3 ) and 

R( x ) = x -^-^F 1 (x 13 ) + x-( 1+ ^F 2 (x* 3 ), (5.1) 

where 

0</3<l + 7 (5.2) 

and 4>(t), Fi(t) and F 2 (t) have the same period u in t, with nowhere zero. Here we have 
N = 2 in Q) and 

01 = 1 + 7-/3, 2 = l + 7 (5.3) 
in (|3+|)-(p72]). Corresponding to ( |3.6| ), we make the induction hypothesis 

V- = x^^+i-^J- (x 13 ) 



where Uj m {i) has period to in t. Then, by (fZ 



where U. m (t) has period u and the entries 7Ty m in U m are obtained from those in C/ lm by 
the formula 



jim/{dj-di) (i^j). 



Then considering (|3.7| ), we have 

p- x p' m = x-^ +7+1 -«(n^/0)(x' 3 )-cT m x-^^ +1 )(n m /0)(x /3 ) 

= x-^^-^W^x 13 ) +x- (CTm+7+1) H / " 2 (x /3 ) (5.4) 

and hence ( |3.7| ) holds with / = 2. 

We note that the upper bound Q5.2p placed on (3 is a restriction on the frequency 
of oscillations of p and R in this example. The same type of condition was imposed 
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in ( [blast ham 1989 , Example 2.4.1) in connection with the method of repeated diago- 
nalization. When /3 > 1 + 7, the asymptotic solution of ( |1.1[ ) requires transformations 
of an entirely different nature from those based on ( |2.1| ) and ( |2.8| ) ( |Eastham| |1992| ). 
( |Sultanaev| |P9l^ ). 

In order to discuss the output of Algorithm 4.1 for this example, we have to choose 
specific values of (3 and 7, so that part (b) can generate the list of values a m . We make 
the simple choice (3 = 7 = 1, so that Q\ = 1 and 9 2 = 2 in ( |5.3|) . Then, by fl3T3|) , a m = m. 
Also, by (|5.1| ) and ( fh4] ), we have 

l Fi(x), 



Vuix) 



x 



V 21 (x) 



X 



W lm {x) = x-^Wi(x), W 2m (x) = x^ m+2 ~>W 2 (x) 

in the notation of section 3. 

We have noted in section 4 that S m+ i depends on Dk, Pk,Wjk and 14 (1 < k < m). 
However, the formulae for the S m+ i can often be simplified by expressing them in terms 
of the earlier Sk- This reduces the number of terms in the formulae with a consequent 
reduction in the computational effort required. We now give the output for S 2 , S3 and S4: 



S 2 
S3 
Si 



V 11 Pi+T 1 + V 12 -Wn 

-PA + V l2 P 1 +T 2 - W 12 - W 2 i 

T 3 - W l3 



We note that, at this stage, these expressions appear no more complex that those com- 
puted in ( [Brown et al.||l99^ ). However, increased difficulties do occur in the evaluation of 
the entries in T3, W\ 2 , W 2 \ and W31 at the next stage when Algorithm 4.2 is implemented. 
We discuss this point further in Example 5.3. The time needed on a Sun SPARC-station 
10 to compute S4 is 1.6 seconds, which is comparable with the comparative time 1.1833 
seconds reported on in ( [Brown et al .| [T996D . The similar times are a reflection of the low 
number of terms that must be manipulated by the symbolic algebra system. As we re- 
marked previously, the output of Algorithm 4.1 at this point consists only of a set of 
symbolic expressions which satisfy non-commutative multiplication. 
The error term 

Ei = -A 3 W 2 3 + A4T4 - A 4 W U A 4 W 2 4 

+ A 1 Pf(V 12 P 1 -W 21 ) 

+ A 2 (-Px (Ti + V 12 - W u ) + V X2 P X - PiVuPi - W 2l ) P 2 

+ A3V31P3 + A4V 41 P 4 

is more complex than that found in ( [Brown et al.||1996|) , which consists of only 11 additive 
terms. This additional complexity in the E 4 term is reflected in the time needed to compute 
norms when the third stage, Algorithm 4.3, is implemented. 
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5.3. Example 3 

Algorithms 4.2 and 4.3 require the input of specific matrices D\ and V\ and, in this 
example, we introduce a special case of Example 5.2 which arises from the n— th order 
differential equation 

y {n \x) - Q(x)y{x) = 0. (5.5) 

Again, Q(x) contains a periodic factor of the form 

Q{x) = X a f(x (3 ) (5.6) 

where f(t) is periodic in t and nowhere zero, with 0</3<l + ^.Asin ( [Brown et al.||1996j , 
section 3), we write ( |5.5| ) in the system form 

Z' = Q 1/n {D + Q'Q- x - x l n C)Z, (5.7) 

where Z has a first component y, D is the diagonal matrix formed by the n— th roots of 
unity, 

D = dg(u u ...,uJn), (5.8) 

and C is constant with 

dgC = -( n - l)(2n)- 1 I. (5.9) 

It follows from (|5.6|) that 

Q'Q-l-l/n = ^-(l+«/n-«(J'J-H/»)(/) + M -(l+«/n)j-l/n^) i 

Hence ( |5.7| ) is the special case of ( |1 . 1| ) and ( |5.1| ) in which 

7 = a/n, = f'\ F x = PfT^C, F 2 = af'^C. 

We now write ( |5.7|) in the form (|2.2|) ( with m = 1), where dgVn = as in fl2.6|) . Thus 
taking the diagonal terms from Fx over to D and using ( ^.9|) , we define 



D x = D-(n- l)(2n)- l /3x~ {l+a/n -^(f'f- l ~ l/n )(x^)I = D + ^(n- l)pl, (5.10) 
where 

and 

Ri = V u + V 21 (=y0, (5.11) 

where 

V xx = -x- a/n np(C - dgC), y 21 = ax- {1+a/n) f- 1/n (x^)C. 



Thus ( |5.10| ) and ( p.ll| ) are our choice of D\ and V\. 



As in the discussion of Example 5.2, we choose the parameter values (3 — 7 = 1, that 
is, a = n and (3 — 1 in (|5.6|). Finally, we must also choose the values of n in order to 
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complete the requirements for implementing Algorithms 4.2 and 4.3. We choose n = 4, in 
which short calculation gives 

( -3 1 + i 1 1 - i \ 

l-i -3 1+i 1 

1 l-i -3 1 + i 

\l+i 1 l-i -3 ) 

as in ( [Brown et al.| |I~996| , Algorithm 6.2). 



c = -l 



The periodic nature of the function introduces additional matrices over the case 
discussed in Example 5.1 and the theory expounded in ( |Brown et al.|[T996|) . As we re- 
marked above, Algorithm 4.1 needs the specific values of (3 and 7 to be available. A 
consequence of this is that all three algorithms must be run for each set of parameter 
values. However the main additional computational difficulties occur in Algorithms 4.2 
and 4.3. In Algorithm 4.2 the extra matrices generated as a consequence of the periodic 
nature of Q must have their entries evaluated, while in Algorithm 4.3 the norm of the 
error matrix, which is considerably more complex than that which occurs in Example 5.1, 
must be evaluated. 

In order to test the performance of the set of algorithms we have chosen to take 

/(x) = 2 + sinx. (5.12) 

However the performance of Algorithms 4.2 and 4.3 is considerably improved if we work 
with a generic function / together with the simplification rule 

f" = 2-f (5.13) 

and the results that we report are based on this latter situation. A further consequence 
of the extra complexity in Q is that, with the CPU power and memory that we have 
available, we can not evaluate the entries in S$. Thus effectively, we can only perform 4 
iterations of Algorithm 1. The time needed on a SPARC 10 workstation to compute the 
entries in S4 is 190 seconds of CPU time. This compares with the 75 CPU seconds that 
was needed in the work reported on in ( |Brown et al.| |1996j ) . The final algorithm, Algorithm 



4.3, deals with the estimation of the sup. norm of the error matrix, in our case E4. This 
involves first applying the Cauchy and triangle inequalities to each matrix component of 
E4 to obtain an upper bound for || E4 || in terms of the norms of its components. An 
estimate for the norm of the inverse matrix A4 is given by ( |4.1| ). 

These norms are estimated by examining and evaluating each component of each matrix. 
In doing this we encounter terms involving p and its derivatives. In order to obtain upper 
bounds for these terms we symbolically compute expressions for them and note that, for 
this example, 

3 > \f(x)\ > 1 (A < x < 00) 

This gives the necessary bound. Again the increase in the complexity of the expressions, 
resulting from the more complex structure of the initial data means that the computation 
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time that is required is increased over ( Brown et al.| 1996 ). It takes some 1110 CPU 
seconds to compute the norm of E± compared with approximately 550 CPU seconds for 
E4 reported in ( [Brown et al.||199^| ). At X = 40 the bound for the norm of this error matrix 
is 1.63099 x 10" 6 . 
We now compute the factor 

I + P {x) 

discussed in ( [4.5| ) and apply this to the asymptotic solution ( |4.4| ) to yield the asymptotic 
solution of ( p.7| ). We mention that in view of the size of the expressions that are generated 
in (|4.5|), we have chosen to evaluate it at x — 40 using 30 digits of accuracy. 



6. Concluding remarks 
6.1. The transformations of Harris and Lutz 

In the introduction, we indicated that the origins of our basic transformation (|1.4|) - (|1.6| ) 
lie in the 1974 paper of Harris and Lutz. In a subsequent paper ( ([Harris fc Lutz| |1977| , 
2.4)), an extension of ( |1.5| ) is also discussed. Whereas ( |1.5| ) can be described as providing 
a first-order approximation to the exact diagonalization of D + R in ( |1.1| ), the extension 
provides a more accurate second-order approximation. These ideas are also discussed in 
(( [Kastham! |I989| ) pp26-8). 

The question arises whether this extension accelerates the process leading to ( |3.9| ), and 
here we indicate why it does not achieve this objective. The essential feature of both ( |1.5| ) 
and the extension in ( |H arris fc Lutz|[T977| , 2.4) is that they are linear algebraic equations 
to determine P. They do not involve P' . Thus, in ( |2.9| ), both the corresponding definitions 
of P m yield a term p~ x P m which, by (|3.7|), contributes expressions Wj m satisfying (^.8|). 
Now, although we have allowed the possibility that W\ m may be zero, there is no reason to 
suppose that it is necessarily zero, and there is therefore nothing to be gained by departing 
from the simplest definition of P m based on (|1 . 5| ) and (1.6). 



6.2. Other computational algorithms 



Here we indicate how our algorithm has a quite different purpose as compared to the 
algorithms of QDora, Crescenzo fc Tournicr| |1982|) and ( Pietrich] |1992|) , and it is convenient 
to refer specifically to the latter. In (|Dietrich| |1992|) , the differential system is 

Y'(x) = x' 1 B(x)Y(x) 

where x can be a complex variable, 

/ 1 \ 



B 



V -b r , 





-62 



-h J 
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and each b in the last row has a Laurent series 

b(x) = (const. )x c (l + aix' 1 + ...) (x — > oo) (6.1) 

with rational c. With oo as an irregular singular point, the solutions of the corresponding 
n — th order differential equation have the asymptotic form 

/(*){1 + (1)} (6.2) 

where the dominant term f(x) comprises the usual logarithmic and exponential factors. 
The algorithm developed by ( Pietrich||1992|) determines f(x) from a knowledge of B. The 



algorithm of ( Delia Dora, Di Crescenzo and Tournier 1982) also allows sub-dominant 
componants of f{x) to be computed. In contrast, our algorithm is concerned with im- 
provements to the o(l) term in (|6.2j ) and the construction of sub-dominant terms as 
explained above in section 1 and at the end of section 4. It is also the purpose of our 
paper to cover classes of coefficients which, unlike (|6.1| ), contain periodic factors as well 
as powers of x. This was the subject of section 5. 
We are grateful to the referees for raising the issues which are covered in this section. 
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